CalcRiskCube <- function(InitialPurchase,YearlyPurchases,YearlyDemand,RecertSchedule,RecertRate) {
  
  # This function is designed to take the input from the PPODS top left graph
  #   and calculate the Likelihood and Cost for the Risk Cube on the right hand side.
  #
  # YearlyPurchases is the purchase input from the top left graph.
  #   - NxY array
  #   - N is the number of parts
  #   - Y is the number of years
  #
  # YearlyDemand is the demand rate times the recert schedule, 
  #   and is the same at all times.
  #   - NxY array
  #==================================================================================
  
  RvR = c(rep(1,3),rep(0.05,22),0.95*!logical(21),rep(0,3),rep(.95,4),rep(0,4),.95,rep(0,14),rep(0.95,7),rep(0,4),.95,rep(0,6),.95,0,0,rep(.95,6),0,.95,0)
  PartCostI = read.csv("PartCost.csv",header=FALSE, sep=",") #, dec=".",,,)  #read in part cost matrix from CSV file
  RepairCost  = PartCostI$V1
  ReplaceCost = PartCostI$V2
  PartCost = PartCostI$V2*RvR + (1-RvR)*PartCostI$V1
#   plot(PartCostI$V2,type="l",col="dark red")
#   lines(PartCostI$V1,type="l",col="blue")
#   lines(PartCost,type="l",col="green")
#   browser()
  
  PLOT = FALSE
  RecertCost  = 69000 #Cost per missile for recertification at Huntsville (provided by program)
  StorageCost = 1000000 #cost per missile per month for missiles to be stored at the facility
  ShelfLife   = 4
  dimInputs = dim(YearlyDemand)
  N       = dimInputs[1]
  OpYears = dimInputs[2]
  
  ExpiredSpares = CalcExpiredSpares(OpYears,N,YearlyPurchases,YearlyDemand,ShelfLife,InitialPurchase)
  if (PLOT) {
    plot(as.numeric(ExpiredSpares[1,]),ylim=c(min(ExpiredSpares),max(ExpiredSpares)))
    title("Expired Spares")
    for (p in 2:N) {
      lines(as.numeric(ExpiredSpares[p,]))
    }
  }
  
  Stock <- CalcStock(OpYears, N, RecertSchedule, YearlyPurchases, YearlyDemand,ShelfLife,ExpiredSpares,InitialPurchase)
  if (PLOT) {
    plot(as.numeric(Stock[1,]),ylim=c(min(Stock),max(Stock)))
    title("Stock")
    for (p in 2:N) {
      lines(as.numeric(Stock[p,]))
    }
  }
  
  Conf <- CalcConf(OpYears, N, Stock, YearlyDemand)
  
  Diff = Stock - YearlyDemand    #Diff is calculated as the difference between the supply of spare parts and the demnad for spare parts.  The need for parts tranlates to missiles that cannot be returned to the fleet
  if (PLOT) {
    plot(as.numeric(Diff[1,]),ylim=c(min(Diff),max(Diff)))
    title("Diff")
    for (p in 2:N) {
      lines(as.numeric(Diff[p,]))
    }
  }
  
  MissileStock <- CalcMissileStock(OpYears, N, Diff, RecertRate, RecertSchedule)
  if (PLOT) {
    plot(as.numeric(MissileStock))
    title("MissileStock")
  }
  
  Cost <- CalcCost(OpYears, N, PartCost, YearlyPurchases, StorageCost, MissileStock, Diff, RecertRate, RecertCost, ReplaceCost,ExpiredSpares,Stock,YearlyDemand)
  
  Risk = c(Cost,Conf)
#   browser()
  return(Risk)
}
#====================================================================================
CalcExpiredSpares <- function(OpYears,N,Purchase,Demand,ShelfLife,InitialPurchase){
  
  ExpiredSpares <- matrix(nrow=N,ncol=OpYears)
  Expired       <- matrix(nrow=N,ncol=OpYears)
  
  for(Part in 1:N){
    for(t in 1:OpYears){
      if(t<4){
        Expired[Part,t]=0
      }
      else if (t==ShelfLife) { # ASK NICK
        Expired[Part,t] = InitialPurchase
                           -Demand[Part,t-3]
                           -Demand[Part,t-2]
                           -Demand[Part,t-1]
                           -Demand[Part,t]
      }
      else if(t==OpYears){
        Expired[Part,t]=0
      }
      else{
        # Expired Parts are what was bought 4 or more years ago minus what was used in those four years.
        Expired[Part,t] = Purchase[Part,t-3]
          -Demand[Part,t-3]
          -Demand[Part,t-2]
          -Demand[Part,t-1]
          -Demand[Part,t]
      }
    }  
  }
  
  ExpiredSpares <- (Expired*(Expired>0))
  return(ExpiredSpares)
}

#====================================================================================

#Stock function to determine how many missiles are in Stock for a given year.  Based on the Demand, purchase schedule, and part shelf life.
CalcStock <- function(OpYears, N, RecertSchedule, Purchase, Demand,ShelfLife,ExpiredSpares,InitialPurchase){
  
  Stock <- matrix(nrow=N,ncol=OpYears)  #initialize the vecotr to contain the number of parts in inventory (Stock)
  
  Stock[,1] = InitialPurchase     #set initial stock value to scale input

  for (t in 1:(OpYears-1)) {

    if (t==1){
      Stock[,t+1] = Stock[,t] +Purchase[,t+1]
    } #ASK NICK
    if (t < ShelfLife){   #run the spares purchase calculation for the first x years 
      #during which time no spares should reach shelf life
      Stock[,t+1] = Stock[,t]-Demand[,t]+Purchase[,t+1]
    }
    if (t >=ShelfLife) {    #run the spares purchase calculation for the number of years in recert vector
      
      Stock[,t+1] = Stock[,t]-Demand[,t]+Purchase[,t+1]-ExpiredSpares[,t]
    }
    Stock[,t] <- Stock[,t]*(Stock[,t]>0) 
  }
  Stock[,OpYears] <- Stock[,OpYears]*(Stock[,OpYears]>0) 
  
  return(Stock)
}
#end Stock

#====================================================================================

#confidence function calculates the "no stock out confidence" based on the calculated demand and current spares inventory
CalcConf<-function(OpYears, N, Stock, Demand) {
  
  PartConf = matrix(nrow=N,ncol=OpYears) #initialize part confidence matrix to store part stock confidence for each part for each year
  Conf = numeric(OpYears)
  Stock4Conf <- Stock    #initialize stock matrix to be used for confidence calculation
  #   DPt=roundUp(Demand,1)
  
  for (t in 1:OpYears) {
    for(Part in 1:N){
      
      Demand[Part,t]->DPt
      Temp= 0 #initialize a temporary vector to house values for confidence equation before performing summation
      
      if (Stock[Part,t] <=0){#add line here to normalize negative stock values to zero.  Purchase cost is included in the underbuy cost (line 236)
        Stock4Conf[Part,t]=0
      }
      if (Stock[Part,t] >50) {
        print("Interesting...")
        browser()
      }
      for (x in 0:Stock4Conf[Part,t]) {
        
        Temp = Temp + (exp(-DPt)*DPt^x)/factorial(x)
      }
      if (Temp == 0) { #error is being caused by negative stock values being fed into the factorial above
        Temp =1
      }
      
      PartConf[Part,t] = Temp
    }
    Conf[t] = min(PartConf[,t])
#     Conf[t] = prod(PartConf[,t])
    
  }
  
#   TotalConf = min(Conf)
  TotalConf = mean(Conf)
  return(TotalConf)
}
#end Conf

#========================================================================================

roundUp <- function(x,to) {    #function used to round up valudes in the Demand and and purchase matricies to integers
  
  y = to*(x%/%to + as.logical(x%%to))
  
  return(y)
}

#========================================================================================

#Missile Stock takes into account insufficient recert rate and insufficient spare parts to calculate the number of missiles that will be unable to be returned to the fleet
CalcMissileStock<-function(OpYears, N, Diff, RecertRate, RecertSchedule) {
  
  MissileStock = numeric(OpYears) #initialize vector to contain the number of missiles in stock for a given year
  UnfixedMissiles = numeric(OpYears)  #initialize vector to contain the number of unfixed missiles due to insufficient spares for a given year
  
  for (t in 1:OpYears) {
    for(Part in 1:N){
      
      UnfixedMissiles[t] = max(abs(Diff)*(Diff<0)) #number of unfixed missiles due to insufficient spares calculated as the max number of parts needed for any part
      
      if (RecertRate>RecertSchedule[t]) {  #find the number of missiles remaining in stock if Throughput and Capacity are sufficient but spares stock is not
        MissileStock[t]<-UnfixedMissiles[t]
      }
      #end insufficient stock stock w/ sufficient throughput & capacity
      else { #throughput and capcity insufficient
        MissileStock[t] = RecertSchedule[t]+UnfixedMissiles[t]-RecertRate
      }
    }
  }
  #end insufficient throughput and capacity
  return(MissileStock)
  
}#end MissileStock

#========================================================================================

#Cost is a function that determines the cost associated with the acquisition of spares in the purchase plan and the recertification costs for the current plan 
CalcCost <- function(OpYears, N, PartCost, Purchase, StorageCost, MissileStock, Diff,RecertRate, RecertCost, ReplaceCost,ExpiredSpares,Stock,Demand){
  
  AnnualOverbuyCost         = matrix(nrow=N,ncol=OpYears)
  AnnualPartCost            = matrix(nrow=N,ncol=OpYears)
  MissilePurchaseCost       = numeric(OpYears)
  AnnualFixedRecertRateCost = numeric(OpYears)
  AnnualCost                = numeric(OpYears)
  MissileCost = 11672592 #cost per missile
  Cost = 0 
  
  MissilePurchaseCost = MissileStock *MissileCost #this is the cost to buy new missiles to keep the availability constant
  
  AnnualFixedRecertRateCost = RecertRate *RecertCost *0.2 #calcuate the cost of a given throughput and capacity as 20% of the recertification cost per missile
  
  for (t in 1:OpYears) {
    
    if(t==OpYears){
      #Cost impact is wasted money used to buy parts that are not needed
#         AnnualOverbuyCost[Part,t] = (ExpiredSpares[Part,t]+Stock[Part,t])*PartCost[Part]
      AnnualOverbuyCost[,t] = (ExpiredSpares[,t] +Stock[,t] -Demand[,t]) *PartCost } else {
      AnnualOverbuyCost[,t] =  ExpiredSpares[,t] *PartCost
    }
    
    AnnualPartCost[,t] = (Purchase[,t] *PartCost) +AnnualOverbuyCost[,t]  #annual cost of initial spares purchase
    
    AnnualCost[t] = sum(AnnualPartCost[,t]) +MissilePurchaseCost[t] +AnnualFixedRecertRateCost #Calculate the annual cost due to under/overbuying spares
  }
  Cost = sum(AnnualCost)
  
#   CostChart = c(sum(AnnualPartCost),sum(MissilePurchaseCost),sum(AnnualFixedRecertRateCost))
#   pie(CostChart,main=t,col=rainbow(length(CostChart)),labels=c("Part","Missile","RecertRate"))
#   
#   plot(colSums(AnnualPartCost,1),type="l",col="dark red")
#   lines(MissilePurchaseCost,type="l",col="blue")
#   lines(AnnualFixedRecertRateCost,type="l",col="green")
#   
#   plot(Purchase[1,],type="l",xlim=c(0,21),ylim=c(0,max(Purchase)))
#   for (p in 2:length(PartCost)) {
#     
#     lines(Purchase[p,],type="l")
#   }
#   filled.contour(1:102, 1:19, Purchase, levels =0:3, xlim=102,ylim=19) #, color.palette=cm.colors,col=color.palette(4))
  
#   browser()
  
#   library(rgl)
#   library(raster)
#   x <- sample(-500:500, 1000, replace = TRUE)
#   y <- sample(-500:500, 1000, replace = TRUE)
#   z <- rnorm(500, 10, 20)
#   df <- cbind(x, y, z)
#   
#   rst <- raster(ncols = 19, nrows = 102, xmn = 0, xmx = 19, ymn = 0, ymx = 102)
#   df <- cbind(Purchase)
#   rst2 <- rasterize(x = df, y = rst, df[, 3])
#   plot3D(rst2)
  
  return(Cost)
}#end Cost function